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ABSTRACT 

Problems on the identification of two dimensional spatial domains arising in the detection 
and characterization of structural flaws in materials are considered. For a thermal diffusion 
system with external boundary input, observations of the temperature on the surface are used 
in a output least square approach. Parameter estimation techniques based on the “method of 
mappings” are discussed and approximation schemes are developed based on a finite element 
Galerkin approach. Theoretical convergence results for computational techniques are given 
and the results are applied to experimental data for the identification of flaws in thermal 
testing of materials. 
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I. INTRODUCTION 


Domain identification problems are important in engineering design and frequently, such 
problems are treated as a branch of the calculus of variations. Discussions usually involve 
nonlinear optimization techniques, optimal control theory, partial differential equations and 
related numerical algorithms. Domain identification for elliptic systems has been studied 
theoretically and numerically by many authors ( for example, see [5], [8], [16], [17], [18] ). 
Most efforts have been concentrated on problems of optimal shape design which are moti- 
vated by numerous applications to structural, airplane, ship design, etc. ( see [17] and the 
references therein ). In this paper, our concern for domain identification is motivated by 
an application that differs somewhat from these shape design problems. Recently, in space 
structures studies, fiber reinforced composite materials have been widely proposed and, as a 
result, demand has grown for assessing the structural integrity of structures made from these 
materials. An important effort on such problems entails non-destructive evaluation methods 
in thermal tomography. These methods involve an attempt to characterize structural flaws 
( e. g., corrosion, cracks, etc. ) which may not be detectable by visual inspection ( see [9] 
for more details ). Mathematically, these problems can be treated as domain identification 
problems. Initial efforts by the authors on the inverse problems of 2-D thermal tomography 
were detailed in [3], [4]. In this paper, we present a more general version of the ideas pro- 
posed in [3], [4] by using the “method of mappings” (see [16], [17] ). Moreover, we report on 
the successful use of the proposed method with both noisy simulated data and experimental 
data from material studies in the Nondestructive Measurement Science Branch at NASA 
Langley Research Center. Other studies closely related to our work have been discussed in 
[9] and [10]. 

To explain our approach, we focus our attention on a 2-D domain identification problem. 
We consider a bounded domain G in two-dimensional Eucledean space as depicted in Fig. 1. 
Let r ( the region of corrosion ) be the section of the boundary where the detection of 
structural flaws is to be attempted via identification of I\ We assume that the boundary of 
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5 1 


Figure 1 . Spatial domain G(q ) bounded by Sr\F(q). 

G is decomposed into two parts dG = S U T, where S is the known surface of the material. 
The system behavior is governed by the 2 -D diffusion equation, 


— - cAu = 0 
at 

in T x G 

(1) 

with the initial and boundary conditions 



13 

II 

O 

"3" 

on G 

( 2 ) 

du 

c w; = f 

on T x S 

(3) 

du 

fa = ° 

on T x r 

(4) 


where c is a given constant representing the thermal diffusivity and where T denotes the 
time interval (t 0} tf) during which the process is observed. In general, the external boundary 
input / can not be measured directly and must also be estimated. For our efforts here we 
parametrize the flux / as a function of an n 0 -dimensional vector q 0 , i. e., 

/ - /(go)- 

We further assume that the corrosion shape at the unknown section T of the boundary is 
specified by an n \ -dimensional parameter vector < 71 , 

r = r(„). 
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The system domain G = G(qi ) is bounded by the partially unknown boundary SUT(qi). In 
treating experimental data, it is often necessary to estimate the initial data from boundary 
observations at time t Q . The resulting estimate will be dependent on q = (qo,qi) through 
both the domain and the input flux. Hence in our presentation the initial state function is 
set as 

u 0 = u Q (q). 


The system output ( i. e. the observations of the temperature ) is assumed to be on a 
subset E of S and, mathematically, the observation is taken as 

y(t,q) = u(t,q)\z. (5) 

The problem treated theoretically and computationally in this paper is that of identifying, 
from output data {y}, the constant parameter vector q determining the geometrical structure 
of the domain G(qi) and the boundary input function /(go)- Let q = (g 0 >9i) be a constant 
parametrization vector among values in a given parameter set Q. Throughout this paper, 
we assume that 

(H-0) The admissible parameter set Q — Qq x Qi is a compact subset of R* 10 x R ni . 

II. DOMAIN IDENTIFICATION BY THE METHOD OF MAPPINGS 

We consider a reference domain C which is a bounded open set in R 2 . The reference 
domain C is taken in the same topology as that of G(gi). Moreover, the domain C is bounded 
by a smooth boundary dC . Then, we introduce the bijective mapping T(gi) from C into 

G(q i), 

x ~ T(q x ) o z 

as depicted in Fig. 2. Mathematically, we make the hypotheses for the mapping T(qi) as 
follows: 
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Figure 2. Transformation mapping from C into G(qi). 

(H-l) The unknown domain G(q\) is given by 

G(?.) = T( qi )(C) 

with existence of f C dC such that 5 = dC/f = S and 

r( 9l ) = T( 9l )(f), 5 = T( 9l )(5) = 5 

where S is independent of 91. 

If we consider a variational formulation similar to that in [ 13 ], the system dynamics ( 1 ) - 
( 4 ) can be described by 

/ A + cVu- V<f>}dx = f f(q 0 )M)dCs for all <j> £ H\G ( gi )) (6) 

JG(qi) at J s 

with 

u(t 0 ) = u 0 (q), ( 7 ) 

where 7 5 denotes the trace operator on S and d(s is a line element on the boundary curve S . 
Using the mapping T{q{], we may obtain the variational formulation defined on the reference 
domain C . 
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Proposition 1: Define 


u(f) = u(t) o T(qi). 

Then the transformed system dynamics can be described by 

du 

)(«(«),*) =!(,)(#) for all <t> G H l (C) (8) 

with 

u(t 0 ) = Uo (q) o T(q 1 ) ( d = Uo(q )) (9) 

where cr(?i )(*, •) and denote, respectively, a sesquilinear form on H l (C) X H l (C) and 

a linear form on H l (C) such that, for <f>, ^ 

= f [cV<t>-(VT(q 1 ))- t (VT(q 1 ))- 1 V'if> (10) 

Jc 

-cV^-{det(VT( 9l ))}“ 1 (VT( gi )r t (Vr( 9a ))- 1 {Vdet(Vr(g 1 ))}^]d 2 

£(?)(<£) d - J s f(<lo){ls<f>)[ys{det(VT(q 1 ))}~ 1 ]d( s (11) 

respectively. 

Proof : The weak formulation of ( 6 ) can be transformed into 

Jc [^ + cv fi . ( vn ft) )-(vr (ft) )-^ det(vr( 9 i))<fa 

= / /(9o)(7s^)^s for j> € ^(C). 

»/S 

Setting (j> = (det(VT(gj))) -1 ^, we have 

jf + jf[cVu^vi'(9,))-‘(vr( 9l ))- 1 w 

- cVu- {det(Vr(„))}' , (Vr(9 1 ))-‘(VT(, 1 ))- , {Vdet(VT(g 1 ))}fl<fr 

= ffMlst) [7s{det(vr( 9l ))} _1 ] ifs 

from which the representation ( 8 ) follows directly. 
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With some tedious calculations, one can readily establish the following useful conditions on 
the sesquilinear form a. 


Theorem 1: Suppose that the mapping T(ij 1 ) satisfies the conditions : 

(C-l) For each q x £ Qi, we have 

T(q ,) € W l „(V) X 
det(VT(„)) £ W'i(CT); 

(C-2) There exists a constant 8 independent of qi such that, for each qi £ Q x , 

det(Vr( gi )) > 8 > 0; 

(C-3) There exists a constant K\ independent of qi such that, for each q\,q\ £ Q i, 


|T(9i) — 7’(9i)lw r i,(C)xWi 1 (C) — ^1 |?i 9i| , 

|det(vr( 9l )) - det (VTCft))!^ (B) < k x | 9x - ^ a | . 

Then, the sesquilinear form a(q x )(-, •) satisfies the following inequalities: There exist positive 
constants a, such that , for <f>, if} £ we have 

— “I'/’Ihmc) - M<f>\h(c) (12) 

l cr (9i)(^) V’)l < ^| < /’Ih 1 (c)|V'Ih 1 (c) (12) 

l <7 (9i)(^)V') _ cr (9i)(^> 1 /')l < /*l9i - 9il \4>\h^{c)\Mh^c)- (14) 

Proof: For convenience in discussions, the explicit form of the mapping x = T(q x )z is 
assumed to be given by 

r*i=w*,«.) ( i5) 

I i 2 = ip 2 {z,qi) 
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Noting that 


we obtain 


VT(„) = 


§£■(«.) &(*) 

L|*-(9i) 


(vr^ortvix ?,))- 1 = {det(vr(„))}- ! 

(I |(9.)) ! + (&(9,)) J -{&(?.)§&(9.)+ S?(9.)^(9.)} 

[ - {^(9.)l^(9.) + §*(9.)S*(9.)} (§£(*»’ + (t-(90) 2 

For economy of notation, we define the functions 

Out*,?') = {(f^(9>)} + {^7< 5,) } 

o 12 (z,9i) = o 2 i(z,9i) 

^dipi d^ 2 sdfa, 

= -l^ (5l W 5l) + ^ (,, W ,l) ) 

022 ( 2 , 91 ) = 


’"W +1^(9.) 


dz ! 


92:1 


b 1 (z,q 1 ) = -a„(z,9 1 )^-{det(VT(g 1 ))}-a 12 ( Z ,9i)^{det(VT(9i))} 
b 2 (z, qi ) = -a 21 (2,9 1 )^-{det(VT(9 a ))}-a 22 (z, gi )^{det(VT(9 1 ))}. 


Then the sesquilinear form ( 10 ) can be rewritten as 

r o % dd>\ d(f > 2 

o(9i)(^i, <t>7.) = y c [c{det(vr(9i))} £ o;_j(z, gi)^--^ 

+c{det(vr ( 9l ))}- 3 £ <n)^rh\dz. 

j<2 ° Z 2 

The principal part of the differential operator becomes 


(16) 


c{det(VT( gi ))}- 2 £ a X3 {z,q x )U } = 

VJ<2 
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This has the lower bound 


C . 1 

- mm < — r, • — 7 — 

2 [a n ( 2 ,gi) < 122(2 


bR 


Under the condition (C-l), there exists a constant R such that 

I dipi ' 


IV'.I < ^ 




< f? for i, j — 1 , 2 . 


Hence we obtain 


c{det(VT ( gi ))}- 2 £ ^ |C | 2 >0 for (6 R 2 , 

*,i<2 


( 18 ) 


(19) 


which means the operator is strongly elliptic. For the second term in ( 16 ), it follows that 


c{det(VT( 9a ))}- 3 i:6 J ( Z , gi )^ 2 

j< 2 U 3 


d(f> 1 


dz x 


+ ^ 2 ( 2 , 91 )! 


<c|det(Vr( 9l ))r 3 |l^(^,9i)l 

Under (C-l), we note that R can be chosen so that we also have 

|det(vr( 9l ))| < R 


d<h 


dz 7 




d_ 

dzi 


{det(VT(gi))} 


< R 


for 


1 , 2 . 


( 20 ) 


From ( 18 ) and ( 20 ), we obtain 

1^(2, 9 i)| < 4ii 3 for j = 1,2. 
Hence, by virtue of (C-2), it follows that 


c{det(VT(< 7 l ))}' 3 E^.9i)l7 i ^ < 

j<2 OZ J 

Consequently, from ( 19 ) and ( 21 ), we can obtain the coercivity property of the sesquilinear 
form. Namely, we obtain 


AcR 3 


6 3 




dzi 


+ 


d<f> 1 


dz 2 


\<t>2\ • 


( 21 ) 


°{qMA) > 


> 


— f V 
4 R 2 Jc §5 

d(f> 

dz t 

— fez 

8 R 2 Jo K f£ 

d<f> 

dzi 


dz — 


4 cR 3 

S 3 


i£ 


d<f> 


dzi 


\<f>\ dz 


+ \4>\ 2 )dz 


c(512i£ 10 + £ 6 ) 

8 ¥r 2 


/ \+f 

Jc 


dz. 
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Hence there exist constants a and A such that ( 12 ) holds. 


To prove the boundedness of c(qi), we note that 


W(qi){<f>u <h)\ < 


/ { det ( VT („))}- 2 E 

Jc t,j<2 aZt 




+c 


/ {dettVTtg^r'EM^)^^ 

•* C j<2 OZ j 


Using ( 18 ) and (C-2), we can estimate the first integral above by 


c / c {det(vr{ 9l ))}- ! E 

A similar bound for the second term follows from the use of the estimate in ( 21 ) and then 
( 13 ) follows. 

To establish the continuity property of ( 14 ), we note that, for any qi,qi G Q i and 
(j>i,(f>2 G ^{C), 

KgiX^i, <t> 2 ) - ^( 9 i)(^i> ^2)1 < 

Q 1 Q 1 

f £ [{det(VT(q l ))}- 2 a, j (z, qt )-{det(VT(q 1 ))}- 1 a, i (z,q 1 )}-^--^-dz 

JC i,j< 2 1 * 

f •£ [{det( VTfft ))} -3 M z > 9t ) - {det(Vr(i,))}- 3 M^9.)] 

JC j< 2 3 

Using ( 18 ) and the condition (C-3), we may argue that 

Mz.gi) - a tJ (z, ?i)| < ARK, \q, - q,\ for i,j = 1,2. 


dz 


< 


2\/2cR 2 


i 9 


I^iIhUC) 1^21 H l (C) 


Similarly, from ( 20 ) and (C-3), we have 

|{det(VT( 9l ))}- 2 - (dettVT^))} 

|{det(VT( 9l ))}- 3 - {det(VT(g 1 ))}" 


2RK, . _ . 

< -j4~m ~ 9il 

3R 2 K, 


< 


£ 6 


ki - I • 
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Hence we have 


IM z >ft) “ M 2 >9 01 < 

Q 

|0ii( 2 ,?i) - a n( z i 9i)| ^{detfVrfai))} 

+ M*,ft)| ^-{det(VT( 9l ))} - ^-{det(VT( 9l ))> 

+ Wi 2 (z,qi) - a i2 (z,qi)\ — {det(VT( 9l ))} 

+ M*.*)I ^{det(VT(„))}- ^-{det(VT($,))} 

< 12R 2 M\ qi - 9l | . 

Thus, we obtain that ( 14 ) follows from 

fa) - °{qi)(<t>u <fe)l 

4 cRK, (3 R 3 V2R 2 3 J\ . . MJI 

< p + —p— + g + Ift - ft I l^il H‘(c) • 

We introduce a Hilbert space 

W(T) = |<6 <*€L : (T;tf‘(C)), ^ 6 i 2 (T; (/f’(C)) - ) 1 

endowed with the norm 

ldU m = jX l#lJn(o) * + /„ % (hi(c)) . *} 

where (J7 1 (C r ))* denotes the dual space of H l (C). Then we have the following result ( see 
[13, Chapter 3, Theorem 1.2] ). 

Lemma 1: In addition to the hypotheses (H-0), (H-l) and the conditions (C-l) - (C-3) in 
Proposition 1, we assume that 

(C-4) For each q £ Q> we have 

u 0 (q) e £ 2 ( c '); 
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(C-5) For each q 0 G Q 0 , we have 


f(qo) € I 2 (5). 


Then, for each fixed q £ Q, the system ( 8 ), ( 9 ) possesses a unique solution u(g) in W(T). 


We next formulate a boundary identification problem. In doing so, we must specify a 
method of data acquisition for the parabolic system ( 8 ), ( 9 ). Let {£p}™ i be given sensor 
locations on the surface S. For each observation point £ p £ S, we introduce the observation 
region EJ, C S such that 

S*=JV©n5 for i — 1,2, • • ■ ,m, 

where N(Q ) denotes a neighborhood of the point £*. Then the observation mechanism for 
the system ( 8 ), ( 9 ) is given by 

y(t,q) = 7i(t)u(t,q) (22) 


where ?f(i) : H l (C) — > R m is given by 

(23) 

for i = 1,2, • • • ,m. Here the hi denote the weighting or mollifying functions for the mea- 
surement equipment ( an infrared camera in the experiments discussed below ). We make 
the following assumption for the observations. 


(C-6) For each i = 1,2, • • • ,m, we have hi 6 L°°(T X EJ,). 

Then recalling ( See [13, Chapter 1, Theorem 3.2] ) that the trace map is continuous from 
H l {C ) to L 2 (S), we have for u(q) € L 2 (T ; fif 1 (C')) that 


y(t,q) = H(t)u(t,q) 


is well defined with y(q) £ L 2 (T ; R m ). We may then define the fit-to-data criterion function 
used in our identification problems by 

J (q) = \ 1 1 IK*.?) -yd(t)\ 2 dt (24) 

Z J to 
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where Jfa(i) denotes the given observed data. Our identification problem may be stated as 
follows: 

(IDP) Find q* E Q which minimizes J(q ) given by ( 24 ) subject to ( 8 ), ( 9 ). 

To insure existence of solutions to (IDP), additional assumptions on the regularity of u 0 
and / are required. We assume: 

(C-7) The mapping q — ► Uq ( q) is continuous from Q to L 2 (C). 

(C- 8 ) The mapping q 0 —* /(go) is continuous from Q 0 to L 2 (S). 

Using continuity properties of the trace map 7 s, lower semicontinuity arguments, and the 
compactness of Q, it is not difficult to prove directly that solutions of (IDP) exist. However, 
these existence results will also follow directly from the approximation considerations in the 
next section ( See Theorem 2 ). 

III. A COMPUTATIONAL METHOD FOR DOMAIN IDENTIFICATION 

In this section, we consider computational techniques for approximate identification prob- 
lems. The approximation scheme we have employed is based on the use of a finite element 
Galerkin approach to construct a sequence of finite dimensional approximating identification 
problems. Let us choose U“=i { as a set of basis functions in H l {C). That is, for each 

N, are linearly independent and \J N 5 pan {^}^ 1 is dense in H^iC). We choose the 

approximation subspaces as 

H n = span{4>i , <f>z, ■ ■ ■ , <f > #}. 

Then, we can define an approximate solution of ( 8 ) by 

^(*. 9 ) = D ™?(t,q)4>? 

i<N 

where the q) are chosen such that, for j — 1 , 2, • • • , iV, 

+ q )j«) = L( ? )(tf) 


12 



with 


Hence, the system ( 8 ), ( 9 ) can be approximated by solving the system 


CV(t) + = F"(q), 

CV(t 0 ) = <(,) 

where 


[c K ]„ 


for i,j — 1,2 ,N 



for i,j = 1,2 ,N 

[F"(»)]i 

= 

for i — 1, 2, • • • , N 

KM)* 

= 

for i = 1,2,--*, TNT 


= <Uo(9),^f> 

for i - 1, 2, • • • , N. 


The corresponding output becomes 




where 'H N (qi) is the m x iV-matrix given by 


[w"(t)] i , = -HiWf 


for i — 1, 2, ■ • • , m; j = 1,2,-**, JV. 


(25) 

(26) 


Then we seek to solve the approximating identification problem given by: 

(AIDP) N Find q N G Q which minimizes 

J N {q) = \y N ( t ,<i)-yd(t)\ 2 dt (27) 

subject to the approximating system ( 25 ), ( 26 ). 

Since it is readily seen that the approximate solution u N (q) ( and hence the maps y N (t, q ) ) 
depend continuously on q, we have that solutions to these approximate problems exist for 
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each N. Our convergence results for the approximate parameter estimates are summarized 
in the following theorem. 

Theorem 2 : Suppose that hypotheses (H-0), (H-l) and (C-l)-(C-8) hold and let q N be a 
solution of the problem (AIDP) N . Then the sequence {g w } admits a convergent subsequence 
with q Nk — > q* as N — > oo. Moreover, q* is a solution of the problem (IDP). 

Proof : Since q N is a solution of the problem (AIDP) N , it is clear that 

J N (q N ) < J N (q) for all q £ Q. 

Thus, if we can argue that for any q M — » q in Q, 

y N (q M )->y(q) in L 2 {T-R m ) asAT,M->oo, 

then, we can obtain the desired inequality 

J(q*) < J{q) for all q £ Q 

by taking limits in ( 27 ). To show this, it suffices to argue that u N (q N ) — ► u(q*) in 
L\T-,H'(C)) as N — * oo for arbitrary sequences {g^} with q N — ► q* as N — ► oo. This 
follows since the output mapping TL is continuous from L 2 {fT\ H^^Cf) to L 2 (T;i2 m ) ( recall 
that 7 5 is continuous from /f 1 (C r ) to L 2 (S ) and hi £ L°°{fT x SJ,) ). 

To argue this convergence of ii N (q N ) to "u(g*) in L 2 (T; /f 1 (C’)), one can use the theory 
developed in [2]. Since Theorem 1 above guarantees conditions (A),(B) and (C) of [2], and 
the density of U nH n in H l {C ) guarantees that (Cl) of [ 2 ] holds, we only need observe that 
( 8 ), ( 9 ) can be written in the form (in H = L 2 {C ) ) 

u(t) = A(q)u(t) + F(q) (28) 

as required in [2] . The only question to address is that of whether our L as given in ( 1 1 ) gives 
rise to an F(q) with values in L 2 (C) ( in the notation of [2], we have if = L 2 (C), V = H l (C) ). 
From ( 11 ) and the continuity of 7 5 : H l (C) — > L 2 (S), we see that <f> — > L(q)(<f>) is in V*. 
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Since V is dense in H = L 2 (C), we can extend ( by continuity ) in the usual manner so that 
the weak system ( 8 ) can be equivalently written in the form ( 28 ). The desired convergence 
results then follow from Theorem 2.3 of [2] along with the representations (2.2) and (2.6) of 
[ 2 ]- 


IV. PRACTICAL IMPLEMENTATION 

In order to implement the identification scheme proposed in Section 3, we need to con- 
struct the mapping operator T(qi). In this section, we treat the special case of identification 
of a boundary shape r(gi) characterized by a simple function. We consider the bounded 
domain G(qi) in R 2 as follows: 

G(qi) = jx E R 2 | 0 < xi < 1, 0 < x 2 < r(i 1 ,g 1 )| 

where x x — > r(x lt qi ) is some parametrized real function which is assumed to characterize the 
unknown part ( the possible corrosive region ) of the boundary r(gi) as depicted in Fig. 3. 
We assume that the boundary of G(qi ) consists of 

dG( qi ) = S\JT( qi ) 

r(9i) = {x|0 < xj < l,x a = r(x i,gi) with r(0,gi) = r(l,g x ) = /} 

U{x|x! = 0,0 < x 2 < i}U{z|Zl = 1,0 <x 2 </} 

S = {x|0 < Xj < 1, x 2 = 0}. 


We use the reference domain C = (0, 1) x (0, l) and we introduce the mapping from C into 

G(qi), 

x = T i < h)( z ) 

(29) 

Then it is easy to verify that this mapping T(g x ) satisfies hypothesis (H-l). Noting that 

1 r '( z l ,91 )*2 

■1 i 


35 i = V’i( z > 9i) = Z 1 

= -4>2(z,qi) = r - ^' ‘ 


VTfoi) 


0 dam} 
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The linear form L(q ) of ( 11 ) in this case has the representation 

l (9)(^)=/ (31) 

Lemma 2: Suppose that 

(C-9) For each q x G Qi? r (9i) belongs to 1^(0, 1); 

(C-10) There exists a constant v\ such that for each gi G Qi, 

r(zi,qi ) > i/i > 0; 

(C-ll) There exists a constant v 2 such that for each qi>q\ G Qi, 

\r(zi,qi) - K 2 i>9i)lwn,(o,i) ^ ^ki - <h\- 

Then the mapping T(<?i) given by ( 29 ) satisfies the conditions (C-l) to (C-3) in Theorem 1. 

From Lemma 2, under the conditions (C-9) to (C-ll), we can apply the domain identifi- 
cation techniques as outlined in the previous sections to the problem posed in this section. 
In the sequel, we consider linear spline approximations of parameterized functions for r and 
/. For Jfe = 0, 1, let L( A Mk ) be the set of piecewise linear splines ( see [6] for details ) with 
the knot sequence A M * = and basis elements {B^ h }. Then we approximate the 

unknown input function and the unknown corrosion shape function by 


/«.») 

Mo —1 

= E 0 + 

Mo-f no— 1 

E qt a °v**°( 0 + 

Mo 

E 0 

(32) 


t=0 

t=M 0 

*=Mo+no 



Afi —1 

Mi +ni — 1 

Mi 



= E o + 

E + 


(33) 


i— 0 i=zMi *=Afi+ru 


respectively. Here we assume that the rji and 8 t are fixed and the g, are to be estimated. In 
order to ensure the hypotheses of Lemma 2, we impose constraints on Q = Qo x Qi with 

Q k = {q k eR nk | q k <q k <<% }, * = 0,1, (34) 
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where qfc and <$ are given lower and upper bounds, respectively. The bounds q f and q^ are 
determined by i/ lt v 2 and upper bounds of the corrosion function r and its derivative. 

In the reminder of this section, we discuss computer implementation of numerical schemes 
for the problem ( AIDP ) N . For the results reported in this paper, we solved the numerical 
optimization problems using a trust region algorithm. Our trust region scheme can be briefly 
stated as follows: Let { q * = (g 0 Jb 9 ifc)}*=i t 2 t - be a sequence generated by this algorithm. At 
the current point q *, we build a “model” of the cost functional ( we usually choose a quadratic 
model ). Then we define a region around where we believe this model to be an adequate 
approximation of the functional. Using this model, we seek a feasible direction so as to 
guarantee a sufficient decrease in the model of cost. Once we obtain the feasible direction, 
the exact cost functional is evaluated at the new point. If its value has decreased sufficiently, 
this new point is accepted and updated as the next iterate, and a new trust region is selected 
for the next iteration. Otherwise the new point is rejected and the current trust region is 
reduced. The advantage of this algorithm is its global convergence properties; namely, this 
algorithm makes it possible to obtain convergence to a critical point (optimal solution), even 
from starting points ( initial guesses ) that are far away from the optimal solution. For 
detailed discussions, we refer to [7] and [14]. 

For the implementation of the trust region algorithm, we used a Fortran software package 
“OPT” created by Dr. Richard Carter of ICASE ( see [7] for details ). Test computations 
were carried out on the Cray-2S at NASA Langley Research Center. 

V. NUMERICAL EXPERIMENTS 

Throughout both our computational and laboratory experiments, we considered corro- 
sion shape identification for steel samples with corrosion like damage ( of varying size and 
corrosion depth ) as depicted in Fig. 4. Since the sample material was steel, the thermal 
diffusivity coefficient c in ( 1 ) was taken as c = 3.18611 X 10 -2 ( (inch) 2 /sec ). The first 

part of this section is devoted to a summary of computational experiments for the purpose of 
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radius = 0. 1 25 



Figure 4. Sample material with corrosion used in one of the experiments 

better understanding our computational method. This is followed by a discussion in which 
the feasibility and validity of our algorithm using laboratory data is demonstrated. 

V.l. Computational Experiments with Simulated Data 

In our numerical experiments, we used a test example constructed as follows: We chose 
a ‘true’ parameter vector q° = ( q q °) Q y generated the corresponding solution numeri- 

cally, added random noise, and then used this as “data” for our inverse algorithm. In this 
procedure, the dimensions of the unknown parameter vectors were taken as dim(qo) — n o = 
5, dim(qi) = rii = 5. The corresponding numbers of knot sequences in ( 32 ) and ( 33 ) 
were chosen as Mq = 4, M 0 = 0 ( i. e., all the coefficients for / were to be estimated ), 
M x = 32, Mi = 14, respectively. The value of the true parameters were chosen as 

q°o = [1.8, 1.8, 1.8, 1.8, 1.8], 

9l ° = [0.0417468, 0.0289693, 0.0250000, 0.0289693, 0.0417468]. 

Figure 5 depicts the true curve approximated by piecewise linear splines. In order to guar- 
antee the conditions (C-9) - (C-ll), the constrained sets Qo and Qi defining Q in ( 34 ) were 
taken as 

Qo = {q 0 eR 5 \ 0.15 < 9 * < 0.25, i= 1,2,3, 4, 5 } 

and 

Qi = {<71 e R s | 0.01 <q[< 0.055, * = 1,2,3, 4,5 }, 
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0.30 0.40 0.50 0.60 0.70 


Figure 5. The corrosion shape and the approximating linear spline curve 
respectively. The initial state function was preassigned as 

u o (q,x) = 0 in G(g x ). 

Then we have 

tio(g) = u 0 (x ) o T(?i) 

= 0 in C. 

To discretize the system model by a ( bilinear spline ) finite element method, the reference 
domain C is divided into a finite number of elements ( K < N ) and a number 

N of npdes defined by {zj = are selected in C. Each element is preassigned 

as an axiparallel rectangle with nodes at the vertices as shown in Fig. 6. The number 
of finite elements and nodes in the numerical experiments reported on here were set as 


Figure 6. The decomposition of finite elements in C 
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K = 128 (= 32 x 4) and N = 165 (= 33 x 5), respectively. The restriction of <f)f to any 
element e* is given by a bilinear polynomial of the form, 




(-S? 


+ “S 1 




for z = (zi, z 2 ) 6 ejt As= 1,2,--*, -RT and i — 1, 2, • • • , JV. 

The coefficients {a^} can be chosen such that each polynomial form satisfies the properties 
of a piecewise bilinear basis function (see e.g., [1, Chapter 5] ). Integration in the element 
matrices C N , A N (qi) and the element vector F N (q ) can be computed numerically by a Gauss- 
Legendre formula. Thus, the state model ( 8 ), ( 9 ) can be solved numerically by an implicit 
scheme with respect to discrete time t = to-\-ih (i = 0, 1, • • ■ , m*) where h = ( tf — The 

initial and final time, and number of time divisions were taken as t 0 = 0 .,</ = 10. and m t = 
80 in our computations. To obtain the output, the weighting functions in the measurement 
operator were given by 

JO. for i — 1,2, ■ •• ,m and 0. < t < 1. 

*(C — | f or { — i f 2, ■ • - , m and 1. < t < 4. 

The number of observation points was set as m = 25 and sensors are located at 
£ = (*i>* 2 )= (^ + 0-125, o) for i = 1,2, •••,25. 


For these test example computations, simulated data were generated by first 

solving the finite element model ( 25 ), ( 26 ) with the same number ( as used in the model 
solution ) of finite elements and nodes. Random noise at various levels from 0% to 50% 
was then added to the numerical solution, thereby producing simulated noisy “data” for the 
algorithm. 

The evaluations of the cost functional its gradient S/ q J N and its Hessian matrix 
are the computationally expensive parts of our algorithm since these involve integration 
of the states w N (t,q) with respect to time t over T. This can be accomplished by using 
the trapezoidal rule in the Newton-Cotes formula. The gradient and the corresponding 
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Hessian are computed using a forward difference scheme and the BFGS secant update with 
safeguarding Hessian approximation ( see [7] ). In each experiment, the optimization routine 
was implemented using rectangular trust regions. 

Table 1 reports the estimated parameter results for the data without noise and with 
5, 10, 15 and 20% relative noise. Figures 7-11 represent the value of cost function, and 
the estimated parameter functions f(z ly q 0 ) and f(z 1 ,g 1 ) which correspond to the estimated 
boundary input and boundary (corrosion) shape, respectively. These and other test compu- 
tations ( see also [4] ) suggest that the methods perform well with reasonable levels of noise 
( 20% or less ). Tests with higher levels of noise ( e. g. 50% ) produced results indicating 
less than satisfactory recovery of reasonable values for the parameters. 

Table 1. True value and estimated values of q 



True Value 

Noise Free 
iter — 61 

5% Noise 
iter =18 

10% Noise 
iter = 17 

15% Noise 
iter =15 

20% Noise 
iter=9 

■ 1 


0.1800 

0.1813 

0.1820 

0.1770 

0.1857 

sps 


0.1800 

0.1787 

0.1778 

0.1809 

0.1738 



0.1799 

0.1832 

0.1847 

0.1859 

0.2038 

HUES 


0.1800 

0.1783 

0.1777 

0.1797 

0.1721 



0.1800 

0.1823 

0.1826 

0.1799 

0.1879 

■9 

0.04175 

0.04165 

0.04403 

0.04479 

0.04444 

0.04156 

99 

0.02897 

0.02896 

0.02814 

0.02885 

0.03195 

0.03869 

9? 

0.02500 

0.02509 

0.02571 

0.02486 

0.02755 

0.03813 

9? 

0.02897 

0.02894 

0.02749 

0.02698 

0.02730 

0.03928 

9 1 

0.04175 

.0.04166 

0.04524 

0.04591 

0.04668 

0.04357 

Ai) 


2.473 

xlO -8 

3.626 

xlO -2 

1.490 

xlO -1 

3.351 

xlO -1 

5.899 

xlO " 1 


( iter = No. of iteration ) 


V.2. Computational Results Using Experimental Data 


In this subsection, we report some of the results of using our estimation procedures with 
experimental data. A schematic of the experimental system is illustrated in Fig. 12. All 
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Figure 7. True function and estimated function (Noise Free) 
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Figure 9. True function and estimated function (10% Noise) 










REMOTE 

SAMPLE 


Figure 12. Block diagram of experimental setup 

experiments were carried out in the laboratories of the Instrument Reseach Division, Nonde- 
structive Measurement Science Branch, NASA Langley Research Center. The experimental 
data consisted of surface temperature distributions for a four second period after a thermal 
source ( heat lamps ) was provided to a sample. In a series of experiments, two types of 
material samples were used : one type of sample fabricated to simulate corrosion similar to 
that shown in Fig. 4 and another type of sample fabricated to simulate no corrosion ( an 
unflawed rectangular sample 0.055 inch by one inch ). The measured experimental data are 
graphed in Fig. 13. 

As readily seen in Fig. 13, the unknown boundary flux / is a time- varying function, 
especially in a neighborhood of the initial time ( this is mainly due to the measurement 
apparutus for the heat source ), whereas, in our algorithm, the boundary input is assumed 
to be a time- invariant function. To handle this discrepancy, noting that the time evolution 
of the experimental data is linearly increasing after an appropriate time period, we set the 
initial time as t 0 = 1.875 sec. Figure 14 depicts the initial time t 0 and the time evolution of 
the experimental data ( for the sample with corrosion ) at the midpoint x = 0.5. 
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(a) with corrosion 



(b) without corrosion 


Figure 13. Experimental data ( surface temperature ) 
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0 

Figure 14. The time evolution of data at x = 0.5 
( sample with corrosion ) 

An important feature of the experimental data as depicted in Fig. 15 is that the magnitude 
of noise is rather high. As demonstrated in the computational experiments of the previous 
section, this could pose difficulties for our estimation procedures. One possible method for 
alleviating this dificculty is to use a data smoothing technique. To this end, we used ICSSCU 
( the IMSL version ) which implements a cubic spline data smoothing. Figure 15 compares 
the smoothed data and the original experimental data at time £ = £ 0 ( = 1.875 ). 

Another important question which arises in using our algorithm is that of how to preassign 
the initial data We only know the surface data 7sUo(= 7s'&o) from the boundary 

surface measurement at time t = t 0 . One possible solution is to approximate the initial state 
by the steady state solution of 

An = 0 in G(qi) 

with the boundary conditions 

^ = o on iXft) 

u — H~ l {to)zd{to) on S , 



.0 1.0 2.0 3.0 4.0 

Time t 
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Figure 15. The smoothed data vs. experimental data at time t - t 0 

( sample with corrosion ) 

where z d (t 0 ) is the observed surface data. We used a finite element method to solve this 
steady state problem for u(qi) and then set Uo(q) = u(qi). 


Inverse problems such as we have considered here may exhibit a lack of continuous de- 
pendence of the estimated parameters on the data. This often leads to serious dificulties in 
computational efforts. It is helpful to use a regularization technique. To this end, for the 
computer implementation of our algorithm using experimental data, we added a regulariza- 
tion term to the cost function, i. e., 


J N (q) = J N (q ) + V 



2 

L*(0,1) 


( see e. g. [11], [15] ). The results of our estimation procedures for one set of experiments 
( the sample with corrosion had dimensions as shown in Fig. 4 ) are given in Table 2 , Fig. 16 
and Fig. 17. Similar results for a second set of experiments for different material samples are 
depicted in Figures 18 and 19 ( in this case, a sample with actual dimensions 0.09375 x 1. 
inch and corrosion radius 0.25 inch was used ). Note that in this case our procedures were 
able to locate a “corrosion” that was not centered in the sample. 
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Table 2. True value and estimated values of q 0 



With corrosion 

Without 

corrosion 

V 

0. 

10. 

0 . 

10. 

q o 

0.2244 

0.2244 

0.1387 

0.1351 

q o 

0.2892 

0.2939 

0.3430 

0.3462 

q o 

0.2267 

0.2316 

0.2451 

0.2521 

q o 

0.2551 

0.2554 

0.2719 

0.2683 

q o 

0.1958 

0.1952 

0.3623 

0.3653 

q } 

0.02718 

0.04388 

0.04466 

0.05500 

9? 

0.03422 

0.03662 

0.05500 

0.05500 

q i 

0.04128 

0.03544 

0.05500 

0.05500 

q t 

0.05135 

0.03778 

0.05500 

0.05500 

q\ 

0.03142 

0.04404 

0.05500 

0.05500 

Ai) 

2.215 

xl0 _1 

2.252 

xl0 _1 

3.222 

xl0 _1 

3.222 

xl0 _1 


VI. CONCLUDING REMARKS 

In the last section we summarized some of our findings using the methods in this paper 
with laboratory data. To date, we have used 7 different samples of material ( all of steel ); 
at least 17 different experiments were carried out and for each of these experiments we 
used our estimation packages with the resulting data. The several results reported in detail 
above are representative of our findings. In all cases, the algorithm performed as well as or 
better than it did in the examples given above. We believe this provides rather conclusive 
evidence that ( i ) structual flaws of this nature can be successfully detected using thermal 
methods, and ( ii ) the mathematical and computational ideas presented in this paper can 
be effectively used in determining the existence, location, and nature of such material flaws. 
^^0 a re currently pursuing both experimental and computational investigations to further 
refine our methods as well as to test other types of materials ( aluminium, etc ) and flaws 
( e. g. delaminations ) with regard to ease and accuracy in detection and characterization 
of flaws using thermal based methods. 
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Cost value 







Iteration 



Figure 18. Estimated results ( sample with corrosion ) 







The theoretical aspects of this paper can be quite properly considered as a nontrivial 
extension of the theory developed in [3], [4]. Here we present a very general theory employing 
the so-called “method of mappings” prominently used in optimal shape design problems ( see 
[16] and [17] ). The theory allows us to treat estimation of unknown boundary segments in 
quite general domains. These domains include as a special case ( as is seen in Section IV ) 
the types of domains treated in [3], [4]. While these more restrictive domains are precisely 
the ones needed to treat the experiments reported on in Section V, the theory of [3], [4] is 
not adequate for some of our current and ongoing investigations of material samples ( with 
delaminations and holes ). Moreover, in our theoretical formulation above we assume that 
the input flux as well as the boundary shape requires identification. This was done because 
we found it impossible to accurately measure the flux as well as the surface temperature in 
experiments with test samples. The theory of [3], [4], which assumes that the flux is given, is 
not adequate to treat the experimental examples discussed in Section V. Finally, the model 
used in [3], [4] includes a boundary transfer coefficient term as well as a radiative loss term in 
the heat equation. We found experimentally that both of these terms appear to be relatively 
unimportant so we have neglected them in the current model and theoretical development. 

The theoretical developments above rely heavily on the framework given in [2], Recently, 
Lamm has, in [12], developed an extension of the results in [2] which could also be used to 
treat the domain identification problems formulated in Sections I and II. The theory in [12] 
involves parameter dependent state spaces H(p) as well as a fixed reference Hilbert space H. 
The example in Section 5 of [12] indicates how the theory developed there could be applied 
to problems involving unknown domains. Briefly, one has a parameterized domain Q p ( anal- 
ogous to our G(qi) above ) and corresponding state spaces H(p) where the norms in general 
also depend on the unknown parameters. The treatment of [12] involves heavy reliance on 
the properties of a “space-changing” map 7 (p>p) • H(p) ► H(p) w ith 7 (p>p) 7 t(p)7e(p) 

where 7e and 7 1 are “extension” and “truncation” maps satisfying certain regularity condi- 
tions. ( When applied to our problems, these are cordinate change maps. ) In [12], the focus 
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is on the relationship between H(p ) and H(p ) through the maps 7 , j t and 7 e as opposed 
to one on the relationship between the parameter dependent spaces and a fixed reference 
space. The developments in [12] could offer some theoretical advantages ( the computational 
procedures would still best be carried out as presented above in Sections III and IV ) in 
that conditions ( 12 ), ( 13 ) and ( 14 ) of Theorem 1 above can be verified directly in the 
variable domain formulation of the problem. In this case one trades the tedium involving the 
coefficients of the transformed system (8 ), ( 9 ) above for some tedium involving changing 
( but equivalent ) norms and verification of regularity properties for the “space-changing” 
maps ( conditions related to (C-l)-(C-3) of our Theorem 1 above ). We refer the interested 
reader to the example of Section 5 of [ 12 ] for details. 
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